library(patchwork)
library(mcmcplots)
library(ggmcmc)
library(tidybayes)
library(tmap)
tmap_mode("plot")
library(terra)
library(sf)
sf::sf_use_s2(FALSE)
library(tidyverse)
# options
options(scipen = 999)# Equal area projection
equalareaCRS <- '+proj=laea +lon_0=-73.125 +lat_0=0 +datum=WGS84 +units=m +no_defs'
Latam <- st_read('data/Latam_vector.shp', quiet = T) %>% st_set_crs(equalareaCRS)
Latam_countries <- sf::st_read('data/Latam_vector_countries.shp', quiet = T) %>% st_set_crs(equalareaCRS)
Latam_no_islands <- bind_rows(list(Latam_countries %>%
filter(type!='Indeterminate' & type!='Dependency' & type!='Lease') %>%
mutate(iso_a2=ifelse(name_en=='France', 'GF', iso_a2)) %>%
mutate(name_en=ifelse(name_en=='France', 'French Guiana', name_en)) %>%
filter(!iso_a2 %in% c('EC','SX', 'NL', 'HT', 'DO', 'CU', 'CW', 'AW','BS',
'TT', 'GD', 'VC', 'BB', 'LC', 'DM', 'AG', 'KN', 'JM')),
Latam_countries %>%
filter(name=='Ecuador') %>%
st_cast('POLYGON', quiet=T) %>%
mutate(area=st_area(.)) %>% arrange(desc(area)) %>%
head(n=1) %>% dplyr::select(-area))) %>% st_union()
Latam.raster <- terra::rast('data/Latam_raster.tif')
# IUCN. (2022). Herpailurus yagouaroundi (spatial data). International Union for Conservation of Nature. https://www.iucnredlist.org/species/9948/50653167
yaguarundi_IUCN <- sf::st_read('data/yaguarundi_IUCN.shp', quiet = T) %>% sf::st_transform(crs=equalareaCRS)
yaguarundi_IUCN.raster <- terra::rasterize(x = vect(yaguarundi_IUCN),
y = Latam.raster,
field = 'presence',
fun = 'length', background=0) %>% mask(., Latam.raster)# Presence-absence data
PA_pre <- readRDS('data/PA_pre.rds') %>%
filter(!is.na(env.elev) &!is.na(env.bio_4) & !is.na(env.npp)) # remove NA's
PA_post <- readRDS('data/PA_post.rds') %>%
filter(!is.na(env.elev) &!is.na(env.bio_4) & !is.na(env.npp)) # remove NA's
# Presence-only data
PO_pre <- readRDS('data/PO_pre.rds') %>%
filter(!is.na(env.elev) &!is.na(env.bio_4) & !is.na(env.npp) & !is.na(acce) & !is.na(count)) # remove NA's
PO_post <- readRDS('data/PO_post.rds') %>%
filter(!is.na(env.elev) &!is.na(env.bio_4) & !is.na(env.npp) & !is.na(acce) & !is.na(count)) # remove NA's
PA_pre_post <- rbind(PA_pre %>% mutate(time=0), PA_post %>% mutate(time=1))
PO_pre_post <- rbind(PO_pre %>% mutate(time=0), PO_post %>% mutate(time=1)) The fitted.model is an object of class
rjags. For details on how it was generated see the code at
code/yaguarundi_IDM_run.Rmd.
fitted.model <- readRDS('data/yaguarundi_model_fit2.rds')
# as.mcmc.rjags converts an rjags Object to an mcmc or mcmc.list Object.
fitted.model.mcmc <- mcmcplots::as.mcmc.rjags(fitted.model)
# labels for the linear predictor `b`
L.fitted.model.b <- plab("b", list(Covariate = c('Intercept', 'env.elev', 'env.bio_7',
'env.bio_15', 'env.npp', sprintf('spline%i', 1:16))))
# tibble object for the linear predictor `b` extracted from the rjags fitted model
fitted.model.ggs.b <- ggmcmc::ggs(fitted.model.mcmc, par_labels = L.fitted.model.b, family="^b\\[")
# diagnostics
# ggmcmc::ggmcmc(fitted.model.ggs.b, file="docs/model_diagnostics.pdf", param_page=3)